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Abstract 

The real radical ideal of a system of polynomials with finitely many 
complex roots is generated by a system of real polynomials having only 
real roots and free of multiplicities. It is a central object in compu¬ 
tational real algebraic geometry and important as a preconditioner 
for numerical solvers. Lasserre and co-workers have shown that the 
real radical ideal of real polynomial systems with finitely many real 
solutions can be determined by a combination of semi-definite pro¬ 
gramming (SDP) and geometric involution techniques. A conjectured 
extension of such methods to positive dimensional polynomial systems 
has been given recently by Ma, Wang and Zhi. 

We show that regularity in the form of the Slater constraint qualifi¬ 
cation (strict feasibility) fails for the resulting SDP feasibility problems. 

Facial reduction is then a popular technique whereby SDP problems 
that fail strict feasibility can be regularized by projecting onto a face 
of the convex cone of semi-definite problems. 

In this paper we introduce a framework for combining facial reduc¬ 
tion with such SDP methods for analyzing 0 and positive dimensional 
real ideals of real polynomial systems. The SDP methods are imple¬ 
mented in MATLAB and our geometric involutive form is implemented 
in Maple. We use two approaches to find a feasible moment matrix. We 
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use an interior point method within the CVX package for MATLAB 
and also the Douglas-Rachford (DR) projection-reflection method. 

Illustrative examples show the advantages of the DR approach for 
some problems over standard interior point methods. We also see the 
advantage of facial reduction both in regularizing the problem and also 
in reducing the dimension of the moment matrices. Problems requiring 
more than one facial reduction are also presented. 


1 Introduction 

In breathrough work Lasserre and collaborators PBi, S3] have shown that the 
real radical ideal of real polynomial systems with finitely many real solutions 
can be determined by a combination of SDP and geometric involution tech¬ 
niques. The real radical ideal of a system of polynomials with finitely many 
complex roots is generated by a system of real polynomials only having real 
roots and free of multiplicities. It is a central object in computational real 
algebraic geometry and important as a preconditioner for numerical solvers. 
A conjectured extension of such methods to positive dimensional polynomial 
systems has been given recently by Ma, Wang and Zhi |28l l27j . 

The above approaches use the method of moments and the Semi-definite 
Programming, SDP formulation. In this paper we see that the Slater con¬ 
straint qualification, strict feasibility, fails for the SDP formulation resulting 
in an ill-posed feasibility problem. Our main contribution is to use facial 
reduction to project the problem onto the minimal face to help regularize 
these computations. Our approach provides tools for working with the ideals 
involved, and gathering data on the open problem above. 

1.1 SDP and Facial Reduction 

The SDP formulation of the moment problem is equivalent to Ending X for 
the linear feasibility system 


AX = b, 16, Si, (1.1) 

where S+ denotes the convex cone of k X k real symmetric positive semi- 
definite matrices, and A : S\ —> M m is a linear transformation. The stan¬ 
dard regularity assumption for 0 is the Slater constraint qualification or 
strict feasibility assumption: 

there exists X with AX = b, X e intS+ . (1.2) 
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We let X >z 0, 0 denote X € S+ , 6 int S+ , respectively. It is well known 

that the Slater condition holds generically, e.g., m- Surprisingly, many 
SDP problems arising from particular applications, and in particular our 
polynomial system applications, are marginally infeasible, i.e., fail to satisfy 
strict feasibility. This means that the feasible set lies in the boundary of the 
cone, and even the slightest perturbation can make the problem infeasible. 
This creates difficulties with the optimality and duality conditions as well as 
with numerical algorithms. To help regularize such SDP problems so that 
strong duality holds, facial reduction was introduced in 1982 by Borwein and 
Wolkowicz mm- However it was only much later that the power of facial 
reduction was exhibited in many applications, e.g., mm®- Developing 
algorithmic implementations of facial reduction that work for large classes 
of SDP problems and the connections with perturbation and convergence 
analysis has recently been achieved in e.g., mmmm- 

A polynomial system of equations can be viewed as a linear (or coeffi¬ 
cient matrix) function of its monomials [251 39]. This linear function yields 
part of the system of linear constraints in the SDP formulation of polyno¬ 
mial systems. The convex cone for polynomials are semi-definite moment 
matrices encoding the real solutions of the polynomial equations and certain 
generalized Macaulay structure possessed by the polynomial systems. Re¬ 
markable advances have been recently made in this area [2511391 El which is 
an intersection between optimization and algebraic geometry. In this article 
we establish a framework for using facial reduction for such systems and 
then solving the systems using the regularized smaller SDP. 


1.2 Prolongation projection methods for involutive bases of 
polynomial systems 

We now look at the details in the semi-definite linear constraint AX = b for 
the polynomial systems. Polynomial systems are remarkable, in that many 
of their constraints are hidden. For example consider the degree two system 

o 

x — x — 1 = 0, xy — y — 1 = 0. 


A single prolongation of this system to degree 3 is found by multiplying 


them by each of the variables x and y: 

x(x 2 — x — 1) = 

x(xy - y - 1) = 

y(x 2 — x — 1) = 

y(xy -y - 1) = 



x 2 y — xy — x 
x 2 y — xy — y 
xy 2 -y 2 -y. 
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Projecting in our paper loosely means eliminating higher degree monomials 
in favour of lower degree ones. In the prolonged system we can project the 
system from degree 3 to degree 2 by eliminating the highest degree term x 2 y 


that occurs in the second and third equations of (1.3): 


x 2 y — xy — x = 0 
x 2 y — xy — y = 0 


xy + x = xy + y. 


(1.4) 


Consequently we obtain the new projected (hidden) constraint x = y. This 
process of uncovering the hidden polynomial constraints by prolongation 
and projection is effected numerically through our geometric involutive form 
algorithm which has been implemented in Maple (381 134] . 

We note that familiar methods for linear systems of equations are Gaus¬ 
sian elimination, GE, for exact solutions and singular value decompositions, 
SVD, for least squares solutions. For polynomial systems, the corresponding 
method in the exact case uses Grobner Bases [5J; while in the approximate 
case we use geometric involutive bases [38j. 


1.3 Facial Reduction and SDP methods applied to real rad¬ 
ical ideals of polynomial systems 

A major motivation for our paper is the success of the work of Lasserre 
et al [25] which gives a new symbolic-numeric approach for computing the 
real radical ideal of zero dimensional polynomial systems using geometric 
involution and SDP techniques. Zero dimensional real polynomial systems 
are systems with real coefficients and finitely many complex and real roots. 
Another major motivation is the important work on this topic in [28, E7j 
which conjectures an extension of [25] to positive dimensional real radical 
ideals. Such ideals have associated real solution components (manifolds) of 
dimension > 1. (See also the paper [36] for examples and many references.) 

The real radical ideal, RRI , of our system P is the set of all polynomials 
with the same zero set as P. To give the reader an informal introduction 
to RRIs and their interpretation, consider the simple case of univariate poly¬ 
nomials with real coefficients, n = 1. In particular, a real univariate poly¬ 
nomial p(x) can be factored in real factors (x — aj) and conjugate complex 
factors (x — ai), (x — cig) so that 

p(x) = Uj(x - aj) dj U k (x - a k ) rh {x - a k ) rk , (1.5) 

where dj and r k are the multiplicities of the roots. The real polynomial ideal 
generated by p(x) is the set of polynomials of the form g(x)p(x) where g(x) 
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is any real polynomial. The RRI of p(x) is generated by the polynomial 

q{x) = II j(x — dj ). (1.6) 

In many applications we are only interested in real roots, and the RRI shown 
here discards all the complex roots. Moreover it also discards multiplicities 
which is important in improving conditioning for polynomial solvers. Many 
general polynomial system solvers, that are capable of determining all so¬ 
lutions explicitly or implicitly, compute all complex and real roots first. In 
particular a generic system of n degree d polynomials in n variables generi- 
cally has d n roots and potentially very few roots. Thus the development of 
methods that avoid the calculation of the complex roots and multiplicities 
is important for efficiency of polynomial system solvers. 


1.4 Outline 


Since we use sophisticated results from diverse areas, in Section[2]we present 
basic ideas and objects through simple examples. We give a preliminary 
introduction to moment matrices and also give a preliminary simple illus¬ 


tration of the power of facial reduction in Section 2.3 


In Section [3] we give a condensed and more formal description of geomet¬ 
ric involutive bases and related algorithms. In Section [4] we discuss moment 
matrices and related algorithms. 

In Section [5] we discuss the methods we used to solve our SDP feasibility 
problems. Since the polynomial problems we consider fail strict feasibility, 
we will use facial reduction to regularize them. However standard primal- 
dual interior point semi-definite programming packages do not deliver the 
accuracy required to guarantee facial reduction. This motivates us to use 
Douglas-Rachford (DR) projection/reflection methods. 

In Section [6] we will discuss our implementation of facial reduction. In 
Section [7] we give numerical experiments. Our concluding remarks are in 
Section [HJ 


2 Basic setup and illustrative examples 

This paper uses sophisticated methods from diverse areas. To help the 
reader, we informally introduce the methods of the paper and illustrate them 
by simple examples. This helps emphasize that the operations underlying 
our approach are reasonably straightforward. 
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2.1 Real polynomial systems 

For background and references to real algebraic geometry and semi-definite 
programming see e.g., 

We consider a (finite) system of £ polynomials P = {pi,...,pf} C 

M[xi,...,x n ] = M[x], where M[x] is the set of all polynomials with real 

T 

coefficients in the n variables x = (x\ X 2 ■ ■ ■ x n ) . We let d = deg(P) 
denote the degree of the polynomial system , i.e., the maximum of the degrees 
of the polynomials pj in P. The solution set or variety of P is 

V K {pi, ...,pi) ={if !C n : pj(x) = 0, VI < j < £}. (2.1) 

This is the real variety of P if K = M and the complex variety of P if IK = C. 
The real ideal generated by P = {p\,... ,pi} C M[x] is: 

(P ) R = ijP l, ■ ■ • ,Pe )m = {fiPi + • • • + fePi ■ fj G M[x], VI < j < £}. (2.2) 

Monomials are denoted by x a := xf 1 ■ ■ ■ xf n , where a G N n , N is the set of 
nonnegative integers, and the degree of x a is |a| := ||a||i = a\ + • • • + a n . 
It is clear that the degree of each monomial |a| < d, the degree of the 
polynomial. Then for appropriate coefficients afc jQ , and for each k, 

we sort by total degree of |a| in nondecreasing order , > 

with components of a sorted in lexicographic order. 

We can rewrite the system of i polynomials, P, as 


P = 


a k} a x a :k = 

\cx\<d 


(2.4) 


Throughout this paper, we use graded reverse lexicographic order, which 
orders first by degree and then by reverse lexicographic order. This order 
respects the Cartan class of variables, which is important in our numerical 
determination geometric features of polynomial systems such as those in 
Definition 13.31 


Definition 2.1 (Coefficient matrix of P, C(P)). Let xbe the column 
vector of monomials x a with 0 < |a| < d sorted as in (2.3). Suppose that 
the coefficients a kt0 in (2.4) are similarly sorted. Then define the coefficient 
matrix of P by C{P) = (afe iQ ). 


The following lemma follows immediately. 


6 




Lemma 2.1. With C(P), x 1 --'^ defined in Definition 2.1, we have 

P = C(P)x { - d \ 


with C{P) E M £xAr ( n - d ) and iV(n,d) := 
mials in x^- d \ 


d + n 
d 


is the number of mono- 


The well-known presentation of polynomial systems as linear functions 
of their monomials and the related coefficient matrix and its kernel and 
rowspace has been exploited in 00E21G21I3H and in the historical work by 
Macaulay m- 


Example 2.1. Consider the system of two univariate polynomials 


P = {x 8 — x 4 — 2, x 8 — 3x 4 + 2} C M[x]. (2.5) 


Here the coefficient matrix is given by C(P) in the equations 


C(P)x < '- 8 l = 


-2 0 0 0 -1 0 0 0 1 
2000 -3 0001 


1 

x 


( 2 . 6 ) 


A familiar computation for many readers is to eliminate the polynomials 
using a Grobner basis calculation: x 8 — x 4 — 2 — (x 8 — 3x 4 + 2) = 2x 4 — 4 or 
equivalently x 4 — 2. The original 8 degree polynomials can be discarded since 
they are consequences o/x 4 — 2. In particular x 8 — x 4 — 2 = x 4 (x 4 — 2) + (x 4 — 
2) = (x 4 + l)(x 4 — 2) so it lies in the ideal generated by x 4 — 2. Similarly 
x 8 — 3x 4 + 2 lies in the ideal generated by x 4 — 2 and can be discarded. All 
polynomials in the ideal generated by P are polynomial multiples of the single 
polynomial 


x 4 — 2. 


(2.7) 


It is easy to see that every system of univariate polynomials is equivalent to 
a single univariate polynomial by applying such simple operations. For sys¬ 
tems of multivariate polynomials, such a minimal object is called a Grobner 
basis. Grobner bases have been intensively studied and usually consist 
of several polynomials. We use the geometric involutive form algorithm dis¬ 
cussed in Section^ to obtain a numerically stable cousin of Grobner bases. 
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2.2 Moment matrices and polynomials 

Moment matrices combined with SDP provide a method to discard the com¬ 
plex roots in polynomial systems with finitely many roots, such as the two 
complex roots of x 4 —2 in Example |2.l| above. Here we focus on the construc¬ 
tion of moment matrices. For theoretical background the reader is directed 
to e.g., [21126]. 

A moment matrix is an infinite real symmetric matrix M = (M Qi/ g) with 
indices corresponding to the indices of the monomials a, /3 £ N n . Here a is 
the index for rows and (3 is the index for columns. Without loss of generality, 
we assume that IW^o = 1. 

Definition 2.2 (Moment matrix). Let u = {u a : a £ N n , |a| < d} £ 
be a vector of indeterminates where the entries are indexed corresponding to 
the exponent vectors of the monomials in n variables of degree at most d. 
The degree d moment matrix of u is a N (n, d) x N (n, d) symmetric matrix 
with rows and columns corresponding to monomials in n variables of degree 
at most d, and defined as 


M d {u) = M{u) = [u a+ f}\ ]am < d - 


Given a multivariate polynomial system P C R[x], with d = deg (P) and 
M £ ^N(n,d)xN(n,d) truncated real symmetric moment matrix. The 


linear constraints imposed by P are, see (2.9) below, 


C{P) M = 0, 


where C(P ) is the coefficient matrix function given in Definition 2.1 


Example 2.2 (Moment matrix for univariate example x = (xi)). The mo¬ 
ment matrix in the univariate (n = 1 ) case is the infinite matrix whose 
(a, /3) entry is u a+ p and a, (3 £ N given by: 


M(u) = 


u 0 

u 1 

U2 

U3 

ua 

Ul 

U2 

U3 

ua 

u 5 

U2 

u 3 

ua 

u 5 

u§ 

U3 

ua 

U 5 

U(i 

u 7 

ua 

u 5 

U6 

U7 

u 8 


U 0 = 1 . 


( 2 . 8 ) 


Note that (2.8) is a Hankel matrix. In Example \2.1\ a degree 8 input system 
was reduced to a degree 4 output polynomial P = {x 4 — 2}. Let us associate 
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u a -H> x a . Then we recover the polynomial equation using the coefficient 
matrix as 


C(P)u { < A) = ( -2 0 0 0 1 ) 


This implies that in terms of the solution x: 


( 1 \ 

ui 

U2 
U3 

V «4 ) 


= 0 . 


C{P)x^ i \x { -- i) ) T = ( -2 0 0 0 1 ) 


X 


X 

2 


2 

X 


X 

3 


3 

X 


X 

.,4 


1 rr .4 , 


= 0. (2.9) 


In the SDP-moment matrix approach we impose uq = 1. We note that the 
association u a ■H- x a extends to the formal correspondence x a x@ u a+ p. 
This allows for the construction of the truncated moment matrix to degree 
d = 4 of the polynomial system as: 


/ 

1 

Ul 

V-2 

U3 

U4 \ 


U\ 

U2 

U 3 

U4 

Uq 


U-2 

U 3 

U4 

Uq 

Uq 


113 

U4 

U5 

Uq 

U 7 

V 

114 

U 5 

Uq 

U 7 

U& ) 


M(u) = 


Appending the linear constraints, we get 

C(P)M(u) = 0 . 


( 2 . 10 ) 


( 2 . 11 ) 


The linear constraints (2.11) are: 


{U 4 — 2 = 0, lt 5 — 2u\ = 0, Uq — 2U2 = 0, It 7 — 2U3 = 0, Uq ~ 2U4 = 0} (2.12) 

which via the correspondence u a <—> x a is equivalent to {x 4 — 2,x 5 — 2x, x 6 — 
2x 2 ,x‘ — 2x 3 ,x s — 2x 4 }. The equivalent SDP problem here is to find a 
maximal rank generic point u = (u a ) where |a| < 2d in the moment matrix 
with 

M(u) P 0, C(P)M(u) = 0. (2.13) 
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By imposing these simple linear constraints we get an explicit simplified 
moment matrix problem in only three variables: 


1 

Ml 

U2 

u 3 

2 

U] 

M2 

m 3 

2 

2mi 

U2 

u-i 

2 

2u\ 

2m 2 

U 3 

2 

2mi 

2u 2 

2m 3 

2 

2u\ 

2U2 

2m 3 

4 


(2.14) 


We note that the substitution of the linear constraints to simplify the problem 
and reduce the number of variables is equivalent to facial reduction; see 
Section [ 6 ] below. This moment matrix problem in ( |2. 14 ) is then sent to an 
SDP solver to approximately find a vector (u\,U 2 , m 3 ) if possible such that M 
is a positive semi-definite matrix with maximum rank. This solver returns 
an approximation which can be recognized for illustrative convenience as 
(ui,U 2 ,U 3 ) = (0, v/2,0) , uo = 1 ,U 4 = 2. Its associated moment matrix and 
moment matrix kernel are: 


M = 


' 1 
0 

V 2 

0 

2 


0 

V 2 

0 

2 

0 


V2 

0 

2 

0 

2^2 


0 

2 

0 

2 \pl 
0 


2 ' 
0 

2 y/2 
0 
4 



f-A 


(~V2\ 


( 0 ^ 

' 


0 


0 


-y/2 



0 

1 

1 

1 

0 



0 


0 


1 





V 0 ) 


V 0 j 



ker M = span K < 


The kernel yields the generating set of three polynomials 

S = {-2 + x 4 ,-V2 + x 2 ,-V2x + x 3 } 

= {(a/2 + x 2 )(-V2 + x 2 ), —a/2 + x 2 ,x(-y/2 + x 2 )}. 


(2.15) 


The factorization in (2.15) allows a trivial Application of the geometric in- 


volutive form algorithm that yields a geometric involutive basis 

{-V2 + x 2 }. 


(2.16) 


The first and third polynomials in (2.15) are a consequence of—\J2 + x 2 by 

Thus we have a basis of the 


our inclusion test, so are discarded, e.g., 

RRI in (2.16). There are efficient eigenvalue methods that can exploit this 
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geometric form to efficiently numerically compute the roots as eigenvalues 
For such solving methods tailored to the real radical and its 
advantages see [25]. The degree 8 system trivially has two real roots given 
by the polynomial in (2.16), i.e., T2 1 / 4 . 


2.3 A class of univariate geometric polynomials 

In this section we experimentally explore the behavior of our facial reduction 
approach (Facial Douglas-Rachford, or abbreviated as FDR) compared to 
a standard SDP solver (Yalrnip SDP, abbreviated as YSDP) which does 
not use facial reduction. In particular we consider the class of univariate 
geometric polynomials which are the partial sums to odd degree d of the 
geometric series: 


Pd{x ) = 1 + x + x 2 + • • • + x d 1 + x d 


where d = 1, 3, 5,.... Then for odd degree d we have 

p d (x) = (x + 1)(1 + x 2 + • • • + x d ~ 3 + x^ 1 ) 


where the even degree factor 1 + x 2 + • • • + x d ~^ + x d ~ l has only complex 
roots. The d roots are x = exp ^, j = 1, • • • , d, and the non-real roots 
appear in complex conjugate pairs. Consequently a generator for the RRI is 

x + 10 

We solved this class of problems for odd degrees d using both the FDf 0 


method with MATLAB R2013b and the YSDP (Yalrnip SDP, R20140605) 
method. We used a laptop (Windows 8.1, Intel Core(TM) i7-4600U CPU 
@2.10GHz 2.70 GHz, 8 GB RAM, 64-bit OS, x64-based processor). 

The running times (in cpu secs) for both methods are given in Figure 
the range of values for the FDR method is clearly better. 


2.3 


3 Geometric involutive bases 


In this section we introduce the basic objects for geometric involutive bases. 
For details and examples see [33 E!- 

Involutivity originates in the geometry of differential equations. See 
Kuranishi [24] for a famous proof of termination of Cartan’s prolongation 
algorithm for nonlinear partial differential equations. A by-product of these 


1 

2 


We denote the generator of the RRI by {pd (a;))u = {x + 1)k- 

The Facial reduction Douglas Rachford method is presented in Section 


5.2.2 


below. 
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d vs YSDP/FDR times for 1 +x+x 2 +...+x d 



Figure 1: Times (cpu secs) for the FDR method versus the YSDP methods 
applied to Pd(x) = 1 + x + ... + x d for odd degrees 1 < d < 69. The blue 
curve (datal on the left) shows YSDP times and the green curve (data2 on 
the right) shows the significantly better FDR times. 


methods has been their implementation for linear homogeneous partial dif¬ 
ferential equations with constant coefficients, and consequently for polyno¬ 
mial algebraic systems. See |3T] for applications and symbolic algorithms 
for polynomial systems. The symbolic-numeric version of a geometric in- 
volutive form was first described and implemented in Wittkopf and Reid 
m- R was applied to approximate symmetries of differential equations in 
[ 8 ] and to polynomial solving in [3?, 135) [38!. See j iq] where it is applied to 
the deflation of multiplicities in multivariate polynomial solving. 

Definition 3.1. Let P be (as usual) a finite subset o/M[x] of degree d. The 
k-th prolongation of system P is D ( P ) = { x a p : 0 < deg (x a p) < d + k,a G 
N n ,p G P}. 

^ ft 

For example D ( P ) for P = {x 2 — x— 1, xy — y— 1} consists of P together 
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with the 4 polynomials in (1.3). 


Definition 3.2. Given a subspace V of J d := R N ( n ’ d ) and t < d, define 
7r £ (V) as the vectors ofV with the components of degree > d — i discarded. 
Given P C R[x] of degree d define 7 r £ (P) := 7 / ker C{P). The k-th proton- 
gation of the kernel is D k (P) := ker C(D P). 

See for example [38j and the published references in [36j for the stable 
numerical implementations of this paper’s operations using SVD methods. 
In Remark 3.5 of (36i we discuss how prolongation and projection can equiv¬ 
alently be computed in the kernel or rowspace, and how polynomial gener¬ 
ators can always be extracted. Underlying this is a 1 to 1 correspondence 
between the relevant vector spaces (not elements). 

Definition 3.3 (Symbol, class and Cartan involution test). Suppose 
P C R[x] of degree d. The symbol matrix S(P) of P is the submatrix ofC(P) 
corresponding to its degree d monomials. Then the class of a monomial x a 
is the least j such that aj 0. 

Suppose that the columns of S(P) are sorted in descending order by 
class and that it is reduced to Gauss echelon form. For k = 1, 2,..., n define 
the quantities ; as the number of pivots in this reduced matrix of class 
k. In a generic system of coordinates the symbol is involutive if 


Y kfid ] = rank S0P) (3.1) 

fc=i 

Suppose Q C M[a:] has degree d! and a ba sis f or ker C(Q) is given by the 
rows of the matrix B. To extract the (3^ in |3.l| ) at projected degree d < d' 
we first numerically project kerC(<5) onto the subspace J d by deleting the 
coordinates in B of degree > d to give a spanning set B for 7v d '~ d Q. Then 
delete the columns in B corresponding to variables of degree < d to obtain 
a matrix Ad corresponding to the orthogonal complement of the degree d 
symbol. Let ^ be the submatrix of B with columns corresponding to 
variables of class < k. In generic coordinates for k = 1... n: 

Pd ] = ( H + d d 1 1 ~ 1 ) - (rank - rank . 


Then the SVD can approximate the ranks in this equation for carrying out 
the Cartan Test (3.1). 
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Definition 3.4 (Involutive System). A system of polynomials P C M[a;] is 
involutive if dim tv DP = dim P and the symbol of P is involutive. 

Definition 3.5. Let P G M[x] with d = deg P and k, i be integers with k > 0 
and 0 < £ < k + d. Then Tr f D k P is projectively involutive if dim , K t D k P = 
dim 7v i+1 D k+1 P and the symbol of7v e D k P is involutive. 

In [H] we prove that a system is projectively involutive if and only if it 
is involutive. In the following algorithm we seek the smallest k such that 
there exists an £ with n i D k P approximately involutive, and generates the 
same ideal as the input system. We choose the system corresponding to the 
largest such £ < k if there are several such values for the given k. 

Input( Q C M[xi,..., x n ]; tolerance e.); 

Set k := 0, d := deg(Q) and P := ker C(Q) m , 

while / 7 ^ 0 do 

Compute D k (P)', initialize set of involutive systems / := {} ; 
for £ from 0 to (d + k) do 
Compute R := 7 t £ D k (P); 
if R involutive then 
| I:=IU{R} 
end 
end 

Remove systems R from /: D d+k ~^R <£. D k (P)', 
k := k + 1 

end 

Output( Return the polynomial generators of the GIF (R) in I of 

lowest degree d = deg R.) 

Algorithm 1: GIF: Geometric involutive form 

The degree of the geometric involutive basis in our method can be lower 
than that given in H3 [23 since Algorithm [l] updates the generators with 
projections. However in the absence of a proof of determination of the real 
radical the larger moment matrices of {28| can capture new members of the 
real radical in situations where our method has already terminated. 

Additional discussion and examples are given in the long version of our 
work [36] . 
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4 Moment matrices &; algorithms 

In this section we outline algorithms for combining geometric involutive 
form and moment matrix methods; see Definition |2.2| Recall that M = 
M(u ) = denotes the moment matrix indexed by a,/3 for rows and 

columns, respectively. And, d = deg(P), M G ]^N(n,d)xN(n,d)^ anc j ^j ie ij near 
constraints imposed by our system of polynomials P C l[i] are given by the 
coefficient times moment matrix multiplication C{P)M = 0. We let (P) K 
denote the associated polynomial ideal and let 

S 

V(P) R = {f€ R[s] + G (Or > 9j € R[x],m G N + }. 

3 =1 

denote the real radical ideal generated by polynomials P over M. A funda¬ 
mental result [5J that is a consequence of the real nullstellensatz is 

VCP)m = {/(*) € R[x] : f{x) = 0, Vx G V R (P)}. 

Input( P = {pi, ...,Pk} C M[xi,... ,x n ]); 

Set Q 0 '■= P, j := 0; 

while r = d do 

d : = dimker GIF(Qj), Qj+i '■= gen(GIF(<3j)); 

Find u* = u(Q j+ 1 ) G R N ( n ’ 2d l: M{u*) P 0,C(Q j+1 )M(u*) = 0; 
r := rank(M(u*)), Q j+2 := gen(kerM(tt*)); 

j '■= j + 2 

end 

Output(Qj + i C M[xi,..., x n ]; Qj+i is in geometric involutive form ; 

— (Qj+ i)r — (-Or-) 

Algorithm 2: GIF M Method 

Algorithm [2] uses the following subroutines described as Algorithms [3] 
and [4] 

Remark 4.1 (Rank-Dim-Involutive Stopping Criterion). A natural 
termination criterion used in Algorithm^ is that the generators stabilize at 
some iteration and the system is involutive: 

gen(GIF(Q )) = g en(ker M(u*)) andQ involutive where u* = u(Q ) (4.1) 

By \25{ (genQier M(Qj + i))) is a sequence of ideals containing y/(P) ■ We 
get an ascending chain of ideals in a Noetherian ring M[xi,..., x n ]. Hence, 
together with the finiteness of the Cartan-Kuranishi geometric involutive 
form algorithm, Algorithm^ terminates. 
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Input( Q C M[a;i,..., x n \. Set d := deg(Q).); 

Construct the moment matrix to degree 2 d.; 

Use SDP methods to numerically solve for a generic point u* = u(Q) 
that maximizes the rank of the moment matrix subject to the 
constraints C(Q) M(u*) = 0.; 

Output( Return M(it*) >: 0 the moment matrix evaluated at this 
generic point.) 

Algorithm 3: M - Moment Matrix 
lnput( GIF(Q) orkerM(u*) where u* =u(Q).); 

Output(Polynonrial generators corresponding to GIF(Q) or kerM(u*)) 

Algorithm 4: gen 


5 Mathematical background for the projection meth¬ 
ods 

In this section we describe the background for the projection methods for 
finding feasible solutions for the moment problems. An important part of 
these methods is building an efficient matrix representation for the linear 
constraints on the moment matrices resulting from the polynomial systems. 


5.1 Linear constraints for multivariate polynomial moment 
matrices 


Recall that we introduced moment matrices informally by a simple example 
in Section |2.2[ see also Definition |2.2[ Let u a := u ai ,....a„ where a G N n 


and the degree of u a is |a| = a\ + ... + a n . Let (a(<d)) be an array of the 
subscripts a of (u a ) with 0 < |a| < d and sorted as in (2.3). 

Consider a truncated moment matrix M{u) = (u a +p) a p &R N(d,n). The 
generalized truncated moment matrix can be represented as follows, where 


(•) yields 

the addition oi 

the subscripts 

for the fj: 




~(fo(u),f 0 (u)) 

(fo(u 

Jl(u)) 

(M u ),Mu)) . 

■ (Mu] 

,fl(u))~ 


(fi(u),fo(u)) 

(h(u 

Jl(u)) 

(Mu), Mu)) . 

■ (Mu] 

,fl(u)) 

M(u) = 

(f2(u),f 0 (u)) 

(Mu] 

,AW) 

(Mu), Mu)) ■ 

• (Mu 

,fl(u)) 


_(fi(u),fo(u)} 

(fi(u 

Ji(u)) 

(. fi(u),Mu)) ■ 

■ (fi(u 

,fl(u))_ 


He re, ( /p, /i,..., _/)) corresponds to the array (u a ) with 0 < |a| < d sorted as 
in (2.3). We denote the i-th element in (u a ) by u l a . Then fi(u) is u l a . 
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In the univariate case the moment matrices have Hankel structure as 
shown in (2.10). In Table[l]we display a truncated bivariate moment matrix 
partitioned into block submatrices having the same degree. Notice that 


MOO 

MlO 

Moi 

«20 

Mil 

M02 

M30 

M21 

Ul2 

M03 

u 10 

M20 

Mil 

M30 

M21 

Ml2 

M40 

M31 

M22 

M13 

«01 

Mil 

M02 

M 2 1 

Ml2 

M03 

M 3 1 

M22 

M13 

Mo 4 

U20 

M30 

M21 

M40 

M31 

M22 

M50 

M41 

M32 

M23 

Mil 

M21 

Ml2 

M31 

M22 

Ml3 

M41 

M32 

M23 

M14 

M 0 2 

Ul2 

M03 

M22 

M13 

Mo 4 

M32 

M23 

M14 

M05 

M30 

M40 

«31 

M50 

M41 

M32 

M60 

M51 

M42 

M33 

U21 

M 3 1 

«22 

M41 

M32 

M23 

M51 

M42 

M33 

«24 

Ul2 

M22 

M13 

M32 

M23 

M14 

M42 

M33 

M24 

Ml 5 

. Mq3 

M13 

Mq 4 

M23 

M14 

Mq 5 

M33 

M24 

M15 

Mq 6 . 


Table 1: A truncated bivariate moment matrix partitioned into block sub¬ 
matrices having the same degree. 

the matrix in Table |T] is not Hankel. However each of its block matrices is 
rectangular Hankel; though even this feature is lost for multivariate moment 
matrices in more than two variables. 

As mentioned above, without loss of generality we assume that uoo = 1- 
As an abbreviation, we may denote M = M(u ) = M r j(u). 

Besides being a symmetric matrix, the moment matrix also has other 
linear constraints among its entries. One can easily see these constraints 
in the truncated univariate matrix (2.10) and bivariate matrix in Table [I] 
An important requirement of our projection methods is to maintain these 
constraints. For example, in the bivariate case above, the matrix elements 
M(u) i 4 = M(u )22 = U 20 are equal. 

We now outline a simple algorithm to find a non-redundant matrix rep¬ 
resentation of these constraints. To list these constraints we start from the 
first row and traverse the matrix from left to right across the rows and then 
traverse the rows from top to bottom. Note also that we only need examine 
entries above the main diagonal since the matrix is symmetric. 


For (2.10) the first linear constraint traversing from the first row down¬ 
wards is M(u) i 4 = M(u) 22 . We denote e* as the i-th unit vector and 
Eij = \ (ef e j + ejei). To impose this constraint, we construct matrix 
At = E'2‘2 — E14, where t represents the index of the linear constraints and 
t = 2 in this case. The constraint is then given by 

(At , M) = trace((£ , 22 - Eu)M) = 0 . 
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Since we always assume = 1, we need to set A\ = E\ \ . Here A t is 

called the matrix representative of the t-th linear constraint. The collection 
of all such matrix representatives for a given moment matrix is called the 
matrix representation of the moment matrix structure. 

Algorithm [ 5 ] below determines all the (non-redundant) matrix represen¬ 
tatives of the linear constraints defining the matrix representation of the 
multivariate moment matrix structure. 

Input(d ; n); 

Initialize array T = and T(i ) is the i-th element of T. 

Initialize n array S = (s ) with the same length as (a(<d)) an d 

S(i ) = [(1, i); o;(<d)(i)] where S(i ), a(<d)(*) is the i-th element of S , 

i a (<d))- 

Let m be the length of T, t = 2 and A\ = E\\. 
for i from 2 to m, do 
for j from i to m, do 

if there exists an s = [( g , h); a] £ S such that T{i ) + T(j) = a 

then 

| ^ = Eij — Egh, t = t + 1 

else 

Adjoin a new element s = [(i,j);a] to S where 
« = T(i) + T(j) 

end 

end 

end 

Output( Return an array of matrix representatives {At] where t E <5, 

£ = { 1 , 2 ,..., 77 } and 77 is the total number of the linear constraints.); 
Algorithm 5: Matrix representation of moment matrix structure 

There are no redundant relations produced by this algorithm so we can 
avoid an overdetermined system. 

In what follows for applications to multivariate polynomial systems of 
degree d in n variables we have 

k~ N(n,d) = (^ d+ d n (5.1) 

Our main problem is the following. 

Problem 5.1 (Main Problem). Let B be a given (k + 1) x m matrix of full 
column rank. Find u G M 2fc+1 so that 

B t M(u ) = 0, trace EuM(u) = 1, M[u ) A 0. 
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We denote PL k+l , space of generalized Hankel matrices. That is these 
matrices have the multivariate structure whose matrix representation is com¬ 
puted by Algorithm [5] It is well known that the special case of Hankel ma¬ 
trices are notoriously ill-conditioned. This means that the cone C\PL k+1 
is thin. i.e., it is close to the boundary of S k+1 , e.g., [20116111]. Therefore, 


solving Problem 5.1 using semi-definite programming techniques results in 


numerical difficulties. 

5.2 Methods of alternating projection and Douglas-Rachford 
projection-reflection 

To apply the methods of alternating projection, MAP or Douglas-Rachford 
reflection-projection, we want to express the main Problem 5.1 as an equiv¬ 
alent problem with moment matrix M = M(u): 


A(M) = b, 


B t M = 0, 


M £ S^ +1 . 


(5.2) 


Here the linear transformation A is obtained from Algorithm [5] The follow¬ 
ing Corollary |5.1| provides the details of the system that we want to solve. We 
first apply facial reduction and get a smaller system. Recall from Algorithm 
[5j we get an array of representing matrix A t s where t £ £, £ = { 1,2,..., rj}. 

Corollary 5.1. Let V be (fc+1) x (k+l—m) and satisfy V T V = I, V T B = 0. 
Let At V T A t V,\/t € £. Let A : S k+1 ~ m —> be defined by 


A{M) := ((traceAtM) ( ) 


t/vte£ 


(5.3) 


Then the main Problem 5.1 with VMV 1 = M{u ) is equivalent to (5.2), 
i.e., to 

A(M) = e 1 , M € S k+1 ~ m , 

and we get M(u) = VMV T . □ 

Let L denote the matrix representation for A in the linear constraints 
There are two projections we use to update the current 
we look at Vc, the linear manifold projection. For the 


in Corollary 5.1 
point p c . First. 

linear system Lp = b = e 1 where L has full row rank, we solve the nearest 
point problem min {^\\p — p c ||i : Lp = 6}, i.e., we find the projection onto 
the linear manifold for the linear constraints. We use L', the Moore-Penrose 
generalized inverse of L. The residual and the update are then 


r c = b — Lp c ] 


= Pc + L ] 


P+ = Pc 


r r . 


(5.4) 
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Second, we project the updated symmetric matrix P + = Vc(Pc) = sHMat(p+) 
onto the semi-definite cone using the Eckart-Young Theorem [XHj, i.e., we 
diagonalize and zero out the negative eigenvalues. Here sHMat = sHvec* = 
sHveW 1 is both the adjoint and the inverse mapping. We denote V S k, the 
positive semi-definite projection and get the new positive semi-definite ap¬ 
proximation V S k (P+). 


5.2.1 Method of alternating projections 

The MAP method is particularly simple, see e.g., the recent book [19] . We 
begin with an initial estimate, e.g., P c = al E J\A mk for a large a > 0. 
We then repeat the projection steps in Items EM till a sufficiently small 
desired tolerance is obtained in the norm of the residual. 

1. Evaluate the residual r c = b — Lp c . Use the residual to evaluate the 
linear projection and obtain the update 

Pl = V c (Pc). 


2. Evaluate the positive semi-definite projection using the Eckart-Young 
Theorem and update the current approximation 

PSDP = V s k(PL)- 


3. Update the cosine value in (5.5). Then update P c 


Psdp- 


The (linear) convergence rate is measured using cosines of angles from three 
consecutive iterates 


cos( 0 ) 


/ trace {{P L - P c )*(Psdp ~ P L )) \ 
V \\ p L-P c\\\\Psdp - Pl)W ) 


(5.5) 


5.2.2 Douglas-Rachford reflection method 

Recall the projections defined above Vc,Ps k ■ We want to find, see (5.2), 
P G Q n S^ +1 , where Q := jp 6 5^ +1 : A(P) = &} . 


We now apply the Douglas-Rachford (DR) projection/reflection method [15]. 
(See also e.g., j3[[9].) 
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Using the QR algorithm applied to B and A, we start with an initial 
estimate 


Po P 0 with B'Pq = 0 and (1,1) component = 1. (5.6) 

Define the reflections IZciPpsd '■ 5+ +1 — > S+ +1 using the corresponding 
projections, i.e., 

llc(P) ■= 2Vc{P) - P , Ppsd(P) := 2 Vpsd(P) - P, VP G M mk . 


Initialization: We set our current estimate P c = Pq to satisfy (5.6). 
We calculate the residual Resc = R — A* s2Mat(P c ), set normres = 
\\Resell denote the reflected residual Resreflc = Resc and reflected 
point IZpsd = Pc- 


• Iterate: We continue iterating from this point while normres > toler, 
our desired tolerance. 


• We use Resrefl to project the current reflected PSD point IZpsd onto 
the linear manifold to get the projected point Pp = Rpsd+A^ Resrefl. 
Then we reflect to get our second reflection point IZp = 2* Pc—PpsD 

• At this time we set our new/current estimate for convergence to be 
P c = Pnew = (Pc + P-c)/ 2 - 

• We now project P c to get Ppsd- We check the residual here for the 
stopping criteria normres = ||Resell = \\R ~ APpsn\\- 

• We now calculate the first reflection point IZpsd = 2 * Ppsd ~ Pc and 
update the reflected residual Resrefl = R — As2vec(Pp5£>). 


The Douglas-Rachford projection/reflection method is simply: 

1. Start at an initial point Pq G 5 (( +1 satisfying (5.6) 

2 . Iterate: Pj+i = \(Pj + PpsD(Pc(Pj)), for all j = 0 , 1 ,.... 


Also the basic theorem on the convergence of the sequence II c(Xk) k , 
[9l Thm 3.3, Page 11], carma.newcastle.edu.au/jon/cycDRinfeas.pdf. so the 
residuals of the projections of the iterates on one of the sets have to be used 
for the stopping criteria. We use the residual after the projection onto the 
SDP cone since finding the residual with respect to the linear manifold is 
inexpensive. 
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To check the linear convergence rates we use the cosine of the angles for 
the vectors of successive iterates, i.e., for three successive iterates P c , 'R-psd^P-c 
and 


trace ((' Rpsp — PcTC^psd — Pc)) 
I \{Ppsd - T^c)\\ II Ppsd ~ Pc\\ 


6 Facial reduction implementation 

Our moment problem is a feasibility problem of the form 

B t M(u) = 0, M(u) y 0, (6.1) 

where B is a given matrix and M(u ) is a linear function of the variables 
u. Constraints on M[u ) are described in Section 5.2, where the problem is 
changed to equality form and then uses facial reduction to get the form 

A(P) = b, Ph 0. (6.2) 

This form includes the first step of facial reduction using the matrix B, 
see Corollary 5.1 and (5.3). Here A(P) = (trace AjP) 6 M m , for specific 
symmetric matrices Aj. 

The projection methods behave poorly when Slater condition fails. We 
therefore attempt to apply further steps of facial reduction and reduce sys¬ 
tem (6.2) until a strictly feasible point exists. We use the following theorem 
of the alternative or characterization of a strictly feasible point; see e.g., [IB]- 

3P,A(P) = b,P^ 0 
Z = A*y hO,b T y = 0 =^> Z = 0. 


(6.3) 


No te th at if a Z A 0 can be found satisfying the left part of the bottom half 
of (6.3) and for the top half P P 0, (P) = b, then 


0 = b T y = (A(P),y} = (P,Z) 


PZ = 0 


range P C null Z. 


Therefore, if the full column rank matrix W satisfies range W = Z, then we 
can facially reduce the problem using the substitution P = WPW T , i.e., we 
can restrict the feasibility problem in (6.2) to the face W ■ W T . 

We can implement the test in (6.3) in several ways. We suppose that A 
is the matrix representation of A, i.e., we let p = s2vec(P) and then we have 

Ap = („4s2Mat)(s2vec(P)) = AP , A*y = s2Mat(A T y). 
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One way would be to first evaluate the orthogonal matrix 
find v so that 


T^irb U 


and 


s2Mat (A T (Uv)) t 0, trace A*(Uv) = (A(I) T U)v = 1. 

Alternatively, we solved 

p* := min \(b T y) 2 
s.t. A*y y 0 

trace A* y = 1 


7 Numerical experiments 

7.1 Examples of Ma, Wang and Zhi [28] 

Ma, Wang and Zhi [2H1 GZ! present an approach using Pommaret Bases 
coupled with moment matrix completion to approximate the real radical 
ideal of a polynomial variety. We applied our approach to [28], Examples 
4.1-4.6]. with the results shown in Table [2] In each case we obtained a 
geometric involutive basis which can be independently verified as a geometric 
involutive basis for the real radical. In |28j Pommaret bases are successfully 
obtained for the real radical for these examples. 

Here are the 6 systems of polynomials corresponding to the examples in 

m- 


{xl + X\X2 - X\X 3 - Xl ~ X2 + X 3 , X\X2 + x\ ~ X 2 X 3 - X\ - X 2 + X 3 , 


X\X 3 + X2X 3 - x 3 - X\ - x 2 + X 3 } ( 7 . 1 a) 

{x\ - X2, X\X2 - £3} ( 7 . 1 b) 

{x\ + xl + xl - 2, x\ + x\-x 3 } ( 7 . 1 c) 

{x\ + X2X 3 - x\, X\X 3 + X1X2 — x 3 , x 2 x 3 + xl + x\ — Xl} ( 7 . 1 d) 


{(xi -x 2 ){xi + x 2 ) 2 (xi +x 2 2 + x 2 ), (xi -x 2 )(xi +X2) 2 {x\ + x%)} (7.1e) 
{(xi - x 2 )(xi + x 2 )(xi + x\ + x 2 ), (xi - x 2 )(xi + x 2 )(xj + x\)} (7.If) 


System (7.1a) for |28L Example 4.1]: Our GIF algorithm [l] with input 
tolerance 10 -10 shows that the system is already in geometric involutive 
form. The corresponding Pommaret basis is given in [28; Example 4.1]. 
The Pommaret basis looks different from the system, but is just a linear 


3 This can be implemented in e.g., CVX using the norm function or absolute value 
function for the objective, i.e., we minimize \b T y\ rather than using the squared term. 
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combination of the system’s polynomials to accomplish the Grobner like 
requirement for its highest terms under the term ordering prescribed in the 
problem. The resulting coefficient matrix of this GIF form, is a full rank 
m = 3, 3 X 10 matrix which is input to the FDR algorithm. Since it has 
rank m = 3, one facial reduction yields a reduced (10 — m) x (10 — m ) = 
7x7 moment matrix. Application of the FDR algorithm using the reduced 
moment matrix, yields convergence in 13 iterations and 0.09 secs, with a 
projected residual error of 10 -14 . These statistics are shown in Table[2j The 
reduction in moment matrix size from 10 x 10 to a 7 x 7 matrix is recorded in 
the rightmost column of the Table by the fraction y. Determination of this 
reduced moment matrix then yields the full 10 x 10 moment matrix of rank 
r = 7. Since the dimension of the kernel for GIF form is d = 7 = r Algorithm 
[^terminates with the input system as its output. It can be checked that the 
ideal generated by this system is real radical. Our facial reduction algorithms 
in Section [6] provide checks for the existence of additional facial reductions. 
They show that there are no additional facial reductions for this problem. 
Syste m (|7.1d[ ) for |28L Example 4.4]: This is very similar to the previous 
system (7.1a). As (28] notes the coordinates for this example are not delta- 
regular, which they and we remedy by a linear change of coordinates. We 
show that the original system is geometrically involutive, which is equivalent 
to the determination of a Pommaret basis by [28]. Just as in the previous 
example, we form a 10 x 10 moment matrix from the GIF form, which is 
transformed by one facial reduction to a 7x7 matrix. There are no additional 
facial reductions, and the full moment matrix and its rank r are determined. 
We find that dimension of the kernel for GIF form is d = 7 = r, so Algorithm 
[2] terminates with the input system as its output. It can be verified the the 
output is a GIF form for the real radical of the ideal. 

System (7.1b) for |28L Example 4.2]: This is quite similar to the sys¬ 
tems (7.1b) and (7.Id). Our methods are similarly efficiently applied to 
this system. Our GIF algorithm first applied one prolongation to the second 
system (7.1b) to yield a degree 3 system. After projectiing from this de¬ 
gree 3 system it shows that the resulting degree 2 system is involutive and 
consists of 3 polynomials. This degree 2 system is geometrically equivalent 
to the Pommaret basis found by [28]. This system is simply the original 
2 polynomials, together with their compatibility condition or S-polynomial 
% 2 (xi — X 2 ) — x\(x\X 2 — x%) = X 1 X 3 — x\. Thus the input system R is re¬ 
placed with 7 tD R with corresponding 3x10 coefficient matrix. The resulting 
10 x 10 moment matrix is facially reduced to a 7 x 7 moment matrix. As 
in the previous examples, no new relations are detected in the kernel of the 
next moment matrix, d = r = 7 and the algorithm terminates. It can be 
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verified that the GIF form is a basis for the real radical ideal of the input 
system. 

Unlike the systems (7.la),(7.lb),(7.Id), the remaining three systems 
(7.1c),(7.1e),(7.1f) of [28] lead to new members in the kernel of their moment 


Example 4.3]: Our initial application of FDR 


matrices. 

System (7.1c) for 
showed slow convergence. However a random linear change of coordinates 
applied to the input system R dramatically improved the convergence. Ap¬ 
plying the GIF algorithm we found that DR is involutive and has a 8 x 20 
coefficient matrix. The dimension of its kernel is d = 12. Facial reduc¬ 
tion then reduces the 20 X 20 moment matrix to a 12 X 12 moment matrix 
which has rank r = 7 ^ d so the algorithm has not terminated. The new 
member of the real radical arising in the moment matrix kernel can be alter¬ 
natively derived by hand by elimination of two of the systems polynomials: 
xl+X 2 +xl — 2— {x\ + X 2 ~X 3 ) = X3 + X3 —2 = (x3 + 2)(x3 — 1). Then noting, 
as explained in [28], that only the root X 3 = 1 leads to real solutions. The 
GIF form of degree 2 of the new system is computed. Its coefficient matrix 
is 5 x 10 and has kernel of dimension d = 5. We note that even with the 
change of coordinates the FDR iteration of this second moment matrix did 
not initially converge until we reduced the required projected residual error 
for production of the first moment matrix to 10 -14 . The second moment 
matrix then was computed quickly and accurately as a 10 x 10 matrix which 
is reduced by one facial reduction to a 5 x 5 matrix. Since the rank of 
the moment matrix is r = 5 = d our algorithm has terminated. It can be 
checked that the output is equivalent to that found by [28] and that the 
resulting GIF form is a basis for the real radical. 

System (7.1e) for [28L Example 4.5]: Direct application of Algorithm [2] 
to (7.le) is relatively inefficient. Instead of this approach we consider an al¬ 
ternative subsystem approach which has the potential to be applied to larger 
systems. Exploiting subsystem structure is a long established approach in 
system solving. 

We apply Algorithm [2] to the subsystem consisting of the first polynomial 
of P\ = {x\ — X 2 )(x\ + X 2 ) 2 (xi + x\ + X 2 ) of (7.le). The GIF form of Pi is 
just Pi, and its coefficient matrix is 1 x 21 matrix with a kernel of dimension 
d = 20. The corresponding moment matrix is 21 x 21, which is reduced to 
a 20 x 20 matrix after one facial reduction. It has rank r = 18 7^ d. So 
the algorithm has not terminated, and new members of the real radical are 
identified from the kernel of the moment matrix. The new system is degree 
5 and has 3 polynomials. Algorithm GIF shows that the first projection 
of this system is involutive and is a single fourth degree polynomial. Its 
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coefficient matrix is 1 x 15 and its kernel has dimension d = 14. The FDR 
algorithm produces a 15 x 15 moment matrix which facially reduced to a 
14 X 14 moment matrix. The rank of the moment matrix is r = 14 = d. 
The algorithm terminates to coefficient errors within 10“ 10 with output as 
a single polynomial which is approximately: 


(xi - x 2 ){x\ + x 2 ){xi + xl + x 2 ) 


(7.2) 


It can be checked that ( |7.2[ ) is a geometric involutive basis for the real radical 
for the ideal generated by Pi. 

Similarly we apply Algorithm p] to the first polynomial of (7. le) which is 


given by P 2 = (x\ — x 2 )(x\ + x 2 ) 2 (xf + x%). The algorithm now terminates 
with output as a single polynomial which is approximately: 


(xi - x 2 ){x 1 + x 2 ) 


(7.3) 


This can be verified to be a geometric involutive basis for the real radical 
for the ideal generated by P 2 . 

Then we consider the system 


(xi - x 2 )(xi + x 2 )(x\ + x\ + x 2 ), - x 2 )(xi + x 2 ) 


(7.4) 


The calculation for (7.1f)for Example 4.6 below yields a geometric involutive 
basis which is approximately 


(®i - xl) 


(7.5) 


It can be independently checked that this is a GIF form for the real radical 


of the ideal of (7. le). 


Syste m (|7.1f| ) for |28l Example 4.6]: This concerns the real solution of 
Q i = ( |7.1 / ) = ( |7.1 / ) subject to the constraints x\ > 1, x 2 > 1. Applying 
Algorithm [2] to Q \ yields a geometric involutive basis which is approximately 
x\ — x 2 . This can be indepdently verified to be a geometric basis for the 
real radical of Q\. The statistics of this reduction are given in the table in 
the row labeled as Ex 4.6 Q\. 

To impose x\ > 1, x 2 > 1 we substitute x\ = x\ + 1, x 2 = x\ + 1 and 
reduce the resulting polynomial Q 2 with Algorithm [2] We obtain x\ — x 2 in 
agreement with |28l Example 4.6]. The statistics of this reduction are given 
in Table [2] in the row labeled as Ex 4.6 Q 2 . 
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FDR 

FDR 

FDR 

GIF-FDR its 

GIF 

Mom Mtx redn 

Syst. 

(n,d,p) 

# its 

secs 

proj res err 

(# FR ) 

tol 

factors s(M)/s(M) 

Ex4.1 

(3,2,3) 

13 

0.09 

10- 14 

1 ( 1 ) 

10" 10 

10 

7 

Ex4.2 

(3,2,2) 

28 

0.01 

10- 14 

1 ( 1 ) 

10- 10 

10 

7 

Ex4.3 

(3,2,2) 

888, 238 

2.3, 0.6 

0 

1 

£■ 

0 

1 

2 ( 2 , 1 ) 

10-10 

20 10 

12 ’ 5 

Ex4.4 

(3,2,3) 

346 

0.53 

10- 14 

1 ( 1 ) 

10" 10 

10 

7 

Ex4.5 Pi 

(2,5,1) 

22314,50 

37.6, 0.3 

10 -12 , 10 -14 

2 ( 2 , 1 ) 

10-!0 

21 15 

20 ’ 14 

Ex4.5 P 2 

(2,5,1) 

957, 1 

4.4, 0.1 

10 -12 , 10 -14 

2 ( 2 , 1 ) 

10-!0 

21 6 

20 » 5 

Ex4.6 Q i 

(2,4,1) 

170, 1 

1.0, 0.09 

10 -12 , 10 -14 

2 ( 2 , 1 ) 

10" 10 

21 6 

15 ’ 5 

Ex4.6 Q 2 

(1,4,1) 

484, 1 

1.4, 0.08 

10 -12 , 10 -14 

2 ( 2 , 1 ) 

10" 10 

15 6 

14 ’ 5 

Cyl2d 

( 2 , 2 , 1 ) 

10 

0.19 

10- 15 

1 ( 1 ) 

10-10 

6 

5 

Cyl3d 

(3,2,2) 

33 

0.77 

10- 14 

1 ( 1 ) 

10" 10 

20 

12 

Cyl4d 

(4,2,3) 

142 

8.45 

10- 14 

1 ( 1 ) 

10" 10 

70 
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Table 2: Statistics for the application of GIF and FDR to polynomial 

systems: n = number of variables, d = maximum polynomial degree, p = the 
number of polynomials; s(M), s(M) sizes of moment matrix M and the faciallly 
reduced matrix M, resp. Ex 4.1-4.6 are the 6 examples in MWZ [ 23 ; Cyl2d-Cyl4d 
are the intersecting cylinder examples. 

7.2 Intersecting higher dimensional cylinders 

Consider the systems of polynomials defining the intersection of n — 1 cylin¬ 
ders in M n 

Cylnd ■= x 2 ! + xl - 1 , xf + x\ - 1 , ■ • ■ ,x\ + x 2 n -l. ( 7 . 6 ) 

Application of the GIF algorithm to the systems Cyl n d for n = 2,3,4 show 
that the systems become geometrically involutive after 0, 2,3 prolongations 
respectively. Table [2] shows the statistics for the subsequent application of 
Algorithm [2] to these systems. The algorithm converges quickly and ac¬ 
curately. Indeed it can be independently determined that the it yields an 
geometric involutive basis for the real radical. 

Further it can be determined that the cylinders form a complete inter¬ 
section and the length of the prolongation to make them involutive, can be 
determined from the symbol of the initial system m • The lower degree 
system, is geometrically formally integrable, and it would be interesting to 
develop methods based on such lower degree systems, to determine, whether 
one can rule out new members in the kernel of the moment matrix of the 
prolonged involutive system from such lower degree systems. 

Finally we mention that recently certain so-called critical point methods 
have been developed for determining witness points [HI [22] on real compo¬ 
nents of real polynomial systems. Indeed the method developed in [H] is 
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successful in finding a point on every component, if the ideal is both real 
radical, and forms a regular sequence. Consequently the systems above, the 
real radical is an important property for such solvers. Such a regular se¬ 
quence can be checked by dimension computation, we only need a formally 
integrable system which has lower degree than the involutive system, this 
leads to a smaller size of moment matrix. Other interesting related results 
are given in [29]. 

7.3 Example of Matlab routine FDR 
Example 7.1. We first use the matrix from ((77]) 

Bx T = [2 0 0 0 -l] . (7.7) 

The moment matrix we get is the exactly the same as that in m Equation 
(37)]: 


1.0000 

-0.0000 

1.4142 

-0.0000 

2.0000 ' 

-0.0000 

1.4142 

-0.0000 

2.0000 

-0.0000 

1.4142 

-0.0000 

2.0000 

-0.0000 

2.8284 

-0.0000 

2.0000 

-0.0000 

2.8284 

-0.0000 

2.0000 

-0.0000 

2.8284 

-0.0000 

4.0000 


The nullity/kernel matrix of P is the same as in J2S, Equation (37)] as well: 


2 

0 

0 

0 

-1 


0.026491 

-0.81147 

-0.09366 

0.57379 

0.052982 


-0.23757' 

-0.090484 

0.83995 

0.063982 

-0.47515 


though it is difficult to see from the last two columns. 

To check whether the matrix B± in (7.7) provides the same nullity as the 
nullity of the matrix P, one can look at the following short MATLAB code 
and see that it is so, i.e., the rank is correct and the spans do not change. 


Bl=[ B’ 

sqrt2 0-100 
0 sqrt2 0-10] 

B1 = 

2.0000 0 
1.4142 


0 


0 - 1.0000 
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0 - 1.0000 

0 0 








0 

1.4142 

0 

- 1.0000 

0 


>> 

B1=B1 ’ 






B1 

= 







2.0000 

1.4142 

0 





0 

0 

1.4142 





0 

- 1.0000 

0 





0 

0 

- 1.0000 





- 1.0000 

0 

0 




>> 

K= [null(P) 

Bl] 





K = 

= 







0.8099 

0.4053 

0.1922 

2.0000 

1.4142 

0 


-0.2574 

0.1542 

0.7593 

0 

0 

1.4142 


-0.4913 

0.6222 

-0.2930 

0 

- 1.0000 

0 


0.1820 

-0.1091 

-0.5369 

0 

0 

- 1.0000 


-0.0575 

-0.6426 

0.1110 

- 1.0000 

0 

0 

>> 

svd(K) 







K» svd(K) 

2.8284 

2.0000 

1.4142 

0.0000 

0.0000 

Following is the output during the MATLAB program. Note the quick 
and accurate convergence; though we have to remember this is a tiny problem. 
It took 118 iterations to get 15 decimals accuracy. The moment matrix P 
has the correct rank. 

Starting with new B value 

using [no^VV 1 ] as initial starting point for P 

time for matrix repres. 0.0468003 

Starting while loop for Douglas-Rachford algorithm 


iter 

cos-vecs 

norm-proj.-resid. 

PSD-proj-; 

10 

0.9938 

0.04919 

6.23e-05 

20 

1 

0.005256 

6.377e-05 

30 

1 

0.0004443 

6.188e-05 

40 

1 

3.282e-05 

6.23e-05 

50 

1 

2.167e-06 

0.000109 

60 

1 

1.27le-07 

6.467e-05 
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70 

1 

6.36e-09 

6.551e-05 

80 

1 

2.341e-10 

6.251e-05 

90 

1 

3.539e-12 

6.349e-05 

100 

1 

1.037e-12 

6.439e-05 

110 

1 

1.324e-13 

6.572e-05 

118 

1 

7.531e-15 

6.404e-05 


time for iterations/while loop is 0.0780005 
max cosine value is 1 

checking feas error in DRalg.m using ***projected*** last iterate Rpsd 
error for norm(B , *P) is 0 

8 Conclusion 

SDP feasibility problems typically involve the intersection of the convex cone 
of semi-definite matrices with a linear manifold. Their importance in ap¬ 
plications has led to the development of many specific algorithms. However 
these feasibility problems are often marginally infeasible, i.e., they do not 
satisfy strict feasibility as is the case for our polynomial applications. Such 
problems are ill-posed and ill-conditioned. 

The main contribution of this paper is to introduce facial reduction, for 
the class of SDP problems arising from analysis and solution of systems 
of real polynomial equations for real solutions. Facial reduction yields an 
equivalent problem for which there are strictly feasible points and which, in 
addition, are smaller. Facial reduction also reduces the size of the moment 
matrices occurring in the application of SDP methods. For example the 
determination of a k x k moment matrix for a problem with m linearly inde¬ 
pendent constraints is reduced to a (k — m ) x (k — m ) moment matrix by one 
facial reduction. We use facial reduction with our MATLAB implementation 
of Douglas-Rachford iteration (our FDR method). In the case of only one 
constraint, say as in the case of univariate polynomials, one might expect 
that the improvement in convergence due to that facial reduction would be 
minor. However we present a class of geometric univariate polynomials of 
odd degree, where one such facial reduction combined with DR iteration, 
yields the real radical much more efficiently than the standard interior point 
method Yalmip. The high accuracy required by facial reduction and also 
the ill-conditioning commonly encountered in numerical polynomial algebra 
m motivated us to implement Douglas-Rachford iteration. 

A fundamental open problem is to generalize the work of [25, [39j to 
positive dimensional ideals. The algorithm of [[281 [27] for a given input real 
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polynomial system P, modulo the successful application of SDP methods at 
each of its steps, computes a Pommaret basis Q: 

yid, 2 (<2)« 2 (P)u (8-1) 

and would provided a solution to this open problem if it is proved that 
(Q)r = y/(P) r- We believe that the work [[28*1 EZ] establishes an important 
feature - involutivity - that will necessarily be a a main condition of any the¬ 
orem and algorithm characterizing the real radical. Involutivity is a natural 
condition, since any solution of the above open problem using SDP, if it es¬ 
tablishes radical ideal membership, will necessarily need (at least implicitly) 
a real radical Grobner basis. Our algorithm, uses geometric involutivity, and 
similarly gives an intermediate ideal, which constitutes another variation on 
this family of conjectures. 

In addition to implementing an algorithm to determine a first facial 
reduction. We also implemented a test for the existence of additional facial 
reductions beyond the first (e.g. in the cases of Examples 4.3 and 4.5 of 
[28]). By using the CVX package or Douglas-Rachford iteration to solve 
for the auxiliary problem, we can determine that if we need a second facial 
reduction by checking whether the optimal value of the auxiliary problem 
is close to 0. So far only moderate improvements in convergence have been 
obtained by our preliminary implementation for construction of additional 
facial reductions. 

Numerical polynomial algebra has been a rapidly expanding and pop¬ 
ular area [40]. It’s problems are typically very demanding, motivating the 
implementation of methods to improve accuracy. For example Bertini, the 
homotopy package developed for numerical polynomial algebra, uses vari¬ 
able precision arithmetic, with particularly demanding problems requiring 
thousands of digits of precision. Consequently this is also a motivation to 
develop higher accuracy methods, such as the FDR method of this paper. 
Manipulations with radical ideals would be a by-product from such work. 

We provided a small set of examples, that illustrate some aspects of 
our algorithms. In Maple all of our examples were executed with Maple’s 
Digits := 15 and the input tolerance := 10” 10 for the GIF algorithm whch 
intensively uses LAPack’s SVD. Accuracy in the projected residual error 
for our tests were between 10 -14 and 10~ 12 . The normalized generators 
obtained for our experiments had coefficients differing less than ICR 10 from 
the exact coefficients. 

Our implementation of auxiliary facial reductions, as still preliminary 
and needs improvement. Even if the real radical is theoretically accessible, 
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the conditioning of the polynomial system, as measured by the sensitivity 
of changes in the solutions to changes in the coefficients, is a significant 
computational affect. So a more detailed study of this aspect is worthwhile. 


32 



References 


[1] A. Alfakih and H. Wolkowicz. Matrix completion problems. In Hand¬ 
book of semidefinite programming, volume 27 of Internat. Ser. Oper. 
Res. Management Sci., pages 533-545. Kluwer Acad. Publ., Boston, 
MA, 2000. |3| 


[2] A.F. Anjos and J.B. Lasserre, editors. Handbook on Semidefinite, Conic 
and Polynomial Optimization. International Series in Operations Re¬ 
search & Management Science. Springer-Verlag, 2011. [6j [8] 


[3] F.J.A. Artacho, J.M. Borwein, and M.K. Tam. Recent results on 
Douglas-Rachford methods. Serdica Mathematical Journal, 39:313-330, 
2013 . m 

[4] S.G. Bartels and D.J. Higharn. The structured sensitivity of 


Vandermonde-like systems. Numer. Math., 62(l):17-33, 1992. 19 


[5] S. Basu, R. Pollack, and M.-F. Roy. Algorithms in Real Algebraic Ge¬ 
ometry, volume 10 of Algorithms and Computation in Math. Springer- 
Verlag, 2 edition, 2006. as 


[6] B. Beckermann. The condition number of real Vandermonde, Krylov 
and positive definite Hankel matrices. Numer. Math., 85(4):553-577, 

2000 . m 


[7] G. Blekherman, P.A. Parrilo, and R.R. Thomas, editors. Semidefinite 
Optimization and Convex Algebraic Geometry. Number 13 in MOS- 
SIAM Series on Optimization. 2012. [3j[6] 

[8] J. Bonasia, F. Lemaire, G.J. Reid, and L. Zhi. Determination of ap¬ 
proximate symmetries of differential equations. Group Theory and Nu¬ 
merical Analysis, 39:249, 2005. 11, 12 [14 


[9] J.M. Borwein and M.K. Tam. A Cyclic Douglas-Rachford Iteration 


Scheme. J. Optim. Theory Appl., 160(1): 1—29, 2014. 20, 21 


[10] J.M. Borwein and H. Wolkowicz. Facial reduction for a cone-convex 
programming problem. J. Austral. Math. Soc. Ser. A, 30(3):369-380, 
1980/81. [| 

[11] J.M. Borwein and H. Wolkowicz. Regularizing the abstract convex 
program. J. Math. Anal. Appl, 83(2):495-530, 1981. [3] 


33 












[12] Y-L. Cheung, D. Drusvyatskiy, N. Krislock, and H. Wolkowicz. Noisy 
sensor network localization: robust facial reduction and the Pareto fron¬ 
tier. Technical report, University of Waterloo, Waterloo, Ontario, 2014. 
in progress. [3] 

[13] Y.-L. Cheung and H. Wolkowicz. Sensitivity analysis of semidefinite 
programs without strong duality. Technical report, University of Wa¬ 
terloo, Waterloo, Ontario, 2014. submitted June 2014. [3 

[14] David A. Cox, John B. Little, and Don O’Shea. Ideals, Varieties, and 
Algorithms. Springer-Verlag, NY, 2nd edition, 1996. 536 pages. [7] 

[15] Jr.J. Douglas and Jr.H.H. Rachford. On the numerical solution of heat 
conduction problems in two and three space variables. Trans. Amer. 
Math. Soc., 82:421-439, 1956. [20] 

[16] D. Drusvyatskiy, G. Li, and H. Wolkowicz. Alternating projections for 
ill-posed semidenite feasibility problems. Technical report, University 
of Waterloo, Waterloo, Ontario, 2014. submitted Sept. 2014. [3] 


22 


[17] M. Diir, B. Jargalsaikhan, and G. Still. The slater condition is generic in 
linear conic programming. Technical report, University of Trier, Trier, 
Germany, 2012. [3] 

[18] C. Eckart and G. Young. A principal axis transformation for non- 


Hermitian matrices. Bull. Amer. Math. Soc., 45:118-121, 1939. 20 


[19] R. Escalante and M. Raydan. Alternating projection methods, volume 8 
of Fundamentals of Algorithms. Society for Industrial and Applied 


Mathematics (SIAM), Philadelphia, PA, 2011. 20 


[20] W. Gautschi and G. Inglese. Lower bounds for the condition number 


of Vandermonde matrices. Numer. Math., 52(3):241-250, 1988. 19 


[21] V.P. Gerdt and Y.A. Blinkov. Involutive bases of polynomial ideals. 


Mathematics and Computers in Simulation, 45(5):519-541, 1998. 12 


[22] Jonathan D Hauenstein. Numerically computing real points on alge¬ 


braic sets. Acta applicandae mathematicae, 125(1): 105—119, 2013. 27 


[23] N. Krislock and H. Wolkowicz. Explicit sensor network localization 
using semidefinite representations and facial reductions. SIAM Journal 
on Optimization, 20(5):2679-2708, 2010. [3] 


34 









[24] M. Kuranishi. On e. cartan’s prolongation theorem of exterior differ¬ 
ential systems. American Journal of Mathematics, pages 1-47, 1957. 

EH 

[25] J.B. Lasserre, M. Laurent, and P. Rostalski. A prolongation-projection 
algorithm for computing the finite real variety of an ideal. Theoretical 
Computer Science , 410(27):2685-2700, 2009. § § [4, 

[26] M. Laurent and P. Rostalski. The approach of moments for polynomial 
equations. In Miguel F. Anjos and Jean B. Lasserre, editors, Handbook 
on semidefinite, conic and polynomial optimization , International Se¬ 
ries in Operations Research & Management Science, 166, pages 25-60. 
Springer, New York, 2012. Em 

[27] Y. Ma. Polynomial Optimization via Low-rank Matrix Completion and 
Semidefinite Programming. PhD thesis, 2012. 0B 

[28] Y. Ma, C. Wang, and L. Zhi. A certificate for semidefinite relax¬ 

ations in computing positive dimensional real varieties. Technical Re¬ 
port arXiv: 1212.4924, KLMM, Academy of Mathematics and Systems 
Science, CAS, 2012. [2} [4| [53| [24| [25| [26j [27| [30j [31] 

[29] Y. Ma and L. Zhi. Computing real solutions of polynomial systems via 
low-rank moment matrix completion. In Proceedings of the 37th In¬ 
ternational Symposium on Symbolic and Algebraic Computation, pages 
249-256. ACM, 2012. [28] 

[30] F.S. Macaulay and P. Roberts. The algebraic theory of modular systems. 
Number 19. University press Cambridge, 1916. [7] 

[31] H.M. Moller and T. Sauer. H-bases for polynomial interpolation and 
system solving. Advances in Computational Mathematics, 12(4) :335— 
362, 2000. 0 [TT] 

[32] B. Mourrain. Isolated points, duality and residues. Journal of Pure 
and Applied Algebra, 117:469-493, 1997. [7] 

[33] B. Mourrain. A new criterion for normal form algorithms. In Applied 
algebra, algebraic algorithms and error-correcting codes, pages 430-442. 
Springer, 1999. 00 

[34] G.J. Reid, P. Lin, and A.D. Wittkopf. Differential elimination- 
completion algorithms for dae and pdae. Studies in Applied Mathe¬ 
matics, 106(1):1—45, 2001. [| 


m 23 30, pT 


11 15 30 


35 














[35] G.J. Reid, J. Tang, and L. Zhi. A complete symbolic-numeric linear 
method for camera pose determination. In Proceedings of the 2003 
international symposium on Symbolic and algebraic computation, pages 


215-223. ACM, 2003. 11 12 


[36] G.J. Reid, F. Wang, and W. Wu. Geometric involutive bases for posi¬ 
tive dimensional polynomial ideals and sdp methods. Technical report, 
Department of Appl. Math., University of Western Ontario, 2014. [4j 

UMEM 

[37] G.J. Reid and L. Zhi. Solving polynomial systems via symbolic-numeric 
reduction to geometric involutive form. Journal of Symbolic Computa¬ 


tion, 44(3):280-291, 2009. 11 12 


[38] R. Scott, G.J. Reid, W. Wu, and L. Zhi. Geometric involutive bases and 
applications to approximate commutative algebra. In Lorenzo Robbiano 
and John Abbott, editors, Approximate Commutative Algebra, pages 
99-124. Springer, 2010. 


[39] F. Sottile. Real solutions to equations from geometry, volume 57 of 
University Lecture Series. American Mathematical Society, Providence, 

RI, 2011. [2}[3j[6j[30 


[40] Hans J. Stetter. Numerical polynomial algebra. Society for Industrial 
and Applied Mathematics (SIAM), Philadelphia, PA, 2004. [7 11 30 


EH 


[41] A.D. Wittkopf and G.J. Reid. Fast differential elimination in c: The cd- 
iffelim environment. Computer Physics Communications, 139(2): 192— 


217, 2001. 12 


[42] H. Wolkowicz, R. Saigal, and L. Vandenberghe, editors. Handbook of 
semidefinite programming. International Series in Operations Research 
& Management Science, 27. Kluwer Academic Publishers, Boston, MA, 
2000. Theory, algorithms, and applications. [6] 


[43] H. Wolkowicz and Q. Zhao. Semidefinite programming relaxations for 
the graph partitioning problem. Discrete Appl. Math., 96/97:461-479, 
1999. Selected for the special Editors’ Choice, Edition 1999. [3] 


[44] W. Wu and G.J. Reid. Finding points on real solution components and 
applications to differential polynomial systems. In Proceedings of the 
38th international symposium on International symposium on symbolic 


and algebraic computation, pages 339-346. ACM, 2013. 27 


36 
















[45] X. Wu and L. Zhi. Determining singular solutions of polynomial systems 
via symbolic-numeric reduction to geometric involutive forms. Journal 
of Symbolic Computation , 47(3):227-238, 2012. 


12 


[46] Q. Zhao, S.E. Karisch, F. Rendl, and H. Wolkowicz. Semidefinite pro¬ 
gramming relaxations for the quadratic assignment problem. J. Comb. 
Optim., 2(1):71-109, 1998. Semidefinite programming and interior- 
point approaches for combinatorial optimization problems (Fields In¬ 
stitute, Toronto, ON, 1996). [3] 


37 



